Activation entropy of electron transfer reactions 
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We report microscopic calculations of free energies and entropies for intramolecular electron trans- 
fer reactions. The calculation algorithm combines the atomistic geometry and charge distribution of 
a molecular solute obtained from quantum calculations with the microscopic polarization response 
of a polar solvent expressed in terms of its polarization structure factors. The procedure is tested 
on a donor-acceptor complex in which ruthenium donor and cobalt acceptor sites are linked by a 
four-proline polypeptide. The reorganization energies and reaction energy gaps are calculated as 
a function of temperature by using structure factors obtained from our analytical procedure and 
from computer simulations. Good agreement between two procedures and with direct computer 
simulations of the reorganization energy is achieved. The microscopic algorithm is compared to 
the dielectric continuum calculations. We found that the strong dependence of the reorganization 
energy on the solvent refractive index predicted by continuum models is not supported by the mi- 
croscopic theory. Also, the reorganization and overall solvation entropies are substantially larger in 
the microscopic theory compared to continuum models. 



I. INTRODUCTION 

Beginning with work of Marcus on electron transfer 
(ET) between ions dissolved in polar solvents 0, the 
understanding of the dynamics and thermodynamics of 
the nuclear polarization coupled to the transferred elec- 
tron has been viewed as a key component of ET theo- 
ries. The concept of polarization fluctuations as a major 
mechanism driving ET has been extended over the sev- 
eral decades of research from simple molecular solvents 
to a diversity of condensed-phase media of varying com- 
plexity. A significant part of the present experimental 
and theoretical effort is directed toward the understand- 
ing of ET in biology, where this process is a key com- 
ponent of energy transport chains 0, 0, 0, Q • Biological 
systems pose a major challenge to theoretical and compu- 
tational chemistry from at least two viewpoints. First, 
the solvent, including bulk and bound water mem- 
branes, and parts of the polar and polarizable matrix of 
the biopolymer, is highly anisotropic and heterogeneous. 
Second, the geometry of what can be separated as a so- 
lute is often very complex, including concave regions of 
molecular scale occupied by the solvent and regions of 
the biopolymer with a significant mobility of polar and 
ionizable residues. 

Dielectric continuum models accommodate the com- 
plex solute shape by numerical algorithms solving the 
Poisson equation with the boundary conditions defined 
by a dielectric cavity pj. The heterogeneous nature of 
the solvent in the vicinity of a redox site can in principle 
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be included by assigning different dielectric constants to 
its heterogeneous parts [8j. Two fundamental problems 
inevitably arise in this algorithm. The first has been 
well recognized over the years of its application and is 
related to the ambiguity of defining the dielectric cav- 
ity for molecular solutes. This problem is often resolved 
by proper parameterization of the radii of atomic and 
molecular groups of the solute. The second problem is 
much less studied. It is related to the fact that collective 
polarization fluctuations of molecular dielectrics possess 
a finite correlation length which may be comparable to 
the length of concave regions of the solute or some other 
characteristic dimensions significant for solvation ther- 
modynamics. The definition of the dielectric constant 
for polar regions of molecular length is very ambiguous 
and, in addition, once the dielectric constant is defined, 
it is not clear if the dielectric response can fully develop 
on the molecular length scale. 

In addition to the problems in implementing the con- 
tinuum formalism for molecular solutes there are some 
fundamental limitations of the continuum approximation 
itself that may limit its applicability to solvation and 
electron transfer thermodynamics. On the basic level, 
the definition of the molecular cavity should be re-done 
for each particular thermodynamic state of the solvent 
[9llT(il] and/or electronic state of the solute This pre- 
cludes the use of continuum theories with a given cavity 
parametrization to describe derivatives of the solvation 
free energy, e.g. entropy and volume of solvation ^^|- I n 
addition, the calculation of the free energy of ET acti- 
vation requires a proper separation of nuclear solvation 
from the overall solvent response. This problem, actively 
studied by formal theories in the past 0, Q, 0, ITq . 
has been recently addressed by computer simulations 
[T7L ITsl ITsj . Computer simulations have indicated that 
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continuum recipes for the separation of nuclear and elec- 
tronic polarization are unreliable, resulting in too strong 
a dependence of the solvent reorganization (free) energy 
on solvent refractive index. All these limitations call 
for an extension of traditional approaches to solvation 
and ET thermodynamics that would include microscopic 
length-scales of solvent polarization. 

Microscopic theories of solvation are not yet sufficiently 
developed to compete efficiently with continuum models 
in application to solvation of biopolymers. Computer 
simulations provide a very detailed picture of the local 
solvation structure, but their application to solvation of 
large solutes requires very lengthy computations and of- 
ten includes approximations that are hard to control. In 
particular, the dielectric response is very slowly converg- 
ing in simulations and is potentially affected by approx- 
imations made to describe the long-range electrostatic 
forces. Several simulation protocols in which polariza- 
tion response is (partially) integrated out by analytical 
techniques have been proposed 20,21]. Integral equation 
theories have been successfully applied to small solutes 
|22|. but examples of their application to solvation and 
reactivity of large solutes are just a few [2j|. The formu- 
lation of the solvation problem in terms of molecular re- 
sponse functions holds significant promise, as it combines 
the molecular length scale of the polarization response 
with a possibility to accommodate an arbitrary shape of 
the solute IHESIplIllIIIIllIJIMlH Ip, H|. 
A recent re-formulation of the Gaussian model [28| for 
solvation in polar solvents [3^, |3(| shows a good agree- 
ment with simulations of model systems and an abil- 
ity to conform with experiment when applied to ET 
in biomolecules |37j | and ch arg e-transfer complexes |3S| 
and to solvation dynamics 39]. Testing the algorithm, 
referred to as the non-local response function theory 
(NRFT), on model systems for which both computer sim- 
ulations and experiment exist is critical for future appli- 
cations to more complex systems. This is the aim of the 
present contribution. 




FIG. 1: Diagram of the polypeptide donor-spacer-acceptor 
(DSA) complex referred to as complex 1 in the text. 

Testing microscopic solvation theories requires com- 
parison to computer simulations on model, yet realis- 
tic, systems. The current experimental database does 
not provide sufficient accuracy to test various approxima- 
tions entering theoretical algorithms. On the other hand, 



computer experiment offers essentially exact (within the 
accuracy of simulation protocols) integration of the same 
Hamiltonian as the one used in the analytical theory. 
Therefore, the present calculations of the ET thermody- 
namics are compared to recent very extensive Molecu- 
lar Dynamics (MD) simulations ,40] of a donor-spacer- 
acceptor (DSA) complex consisting of transition-metal 
donor (D) and acceptor (A) sites linked by a polyprolinc 
peptide spacer (S) (Fig. [TJ: 

(bpy) 2 Ru 2+ (bpy') - (pro) 4 - 0-Co 3+ (NH 3 ) 5 , 

where in the donor bpy=2,2'-bipyridine and bpy'=4'- 
methyl-2,2'-bipyridyl. The spacer is a polyproline chain 
whose first member (the N-terminus) is connected to the 
bpy' carbonyl, and whose fourth member is terminated 
by a carboxylate moiety bound to the — Co 3+ (NH3)5 ac- 
ceptor. This system, modeling ET in redox proteins, is 
a representative member of a homologous series of DSA 
complexes for which ET rates as a function of tempera- 
ture have been reported [41|. This complex will be re- 
ferred to as complex 1 in the text. 

The analytical NRFT model is shown to agree excep- 
tionally well with MD simulations for complex 1 (Figure 
^) in TIP3P water. In order to provide a rigorous com- 
parison between simulations and analytical theory, the 
set of solute charges employed in the simulations was 
also used in the analytical calculations. In addition, the 
polarization structure factors of TIP3P water were ob- 
tained from separate MD simulations to be used as input 
in the analytical theory. Once the accuracy and robust- 
ness of the analytical procedures are tested on MD sim- 
ulations, the next step is to see if the model is capable 
of reflecting the behavior of real systems. To this end, 
we have developed a parameterization scheme for polar- 
ization structure factors applicable to polarizable polar 
solvents. Once this is done, the theory can be extended 
to calculations at varying thermodynamic conditions of 
the solvent (e.g., temperature) and should generate a set 
of predictions which can be tested experimentally. 

We use the polypeptide DSA to focus on two prob- 
lematic areas of dielectric continuum models: depen- 
dence of the reorganization energy on the solvent po- 
larizability [171 |19| and the entropy of nuclear solva- 
tion [T2I |29|. For both areas there is a fundamen- 
tal, both quantitative and qualitative, disagreement be- 
tween microscopic models and continuum calculations. 
Unfortunately, no experimental evidence on the depen- 
dence of the solvent reorganization energy on solvent re- 
fractive index is available in the literature. There is, 
on the other hand, a limited number of experimental 
[p, El El 13 EE EE E3, EB, E3, El|j and simula- 
tion [52] studies on the entropy of reorganization. Most 
of the available experimental (laboratory and simulation) 
evidence points to a positive reorganization entropy (i.e., 
a negative slope of the reorganization energy vs tempera- 
ture) in polar solvents , in agreement with the prediction 
of microscopic theory [29] and in disagreement with neg- 
ative entropies from continuum calculations [T2II53I] . We 
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are aware, however, of a few measurements performed on 
charged donor-acceptor complexes indicating either zero 
or negative reorganization entropies |42l lisl Efij . Our 
current calculations on complex 1 (Fig. Q give absolute 
values of the reorganization entropy much higher than 
continuum calculations. This great discrepancy calls for 
additional tests of the theory against experimental data, 
which will be a subject of future work. 



II. GOLDEN RULE RATE CONSTANT 



The Golden Rule rate constant of ET is |H3 | 



VET 



^^FCWD(O), 



(1) 



where FCWD stands for the density-of-states weighted 
Franck-Condon (FC) factor 



FCWD(w) 



dt 
2^ 



e iH 2 t/n e -iH 1 t/n\ e -iui 



(2) 



Here, (...)„ is an ensemble average over the nuclear de- 
grees of freedom of the system (denoted by subscript 
"n" ) , which include the manifold of N normal vibrational 
modes of the donor-acceptor complex Q = {qi, . . . qjy} 
and the nuclear component of the dipolar polarization of 
the solvent P n . The ensemble average is carried out over 
the configurations in equilibrium with the initial state. 
Further, Hi (i = 1,2) are the diagonal matrix elements 
of the unperturbed system Hamiltonian H taken on the 
two-state electronic basis {*i,* 2 }: H = (^i\H\^i) 
(i = 1 and i = 2 stand for the initial and hnal electronic 
states, respectively). The sum of H and the perturbation 
V makes the whole system Hamiltonian, H' = H + V , 
and V12 = (vE'i|^| v I'2) is the off-diagonal matrix element 
in the Golden Rule expression. 

The system Hamiltonian of a donor-acceptor complex 
in a condensed-phase solvent can be separated into the 
gas-phase component, H g , the solute-solvent interaction, 
Hq s ("0" stands for the solute, "s" stands for the solvent), 
and the bath Hamiltonian, Hb, describing thermal fluc- 
tuations of the solvent: 



H = H a 



H 



0s 



(3) 



The gas-phase Hamiltonian is the sum of the kinetic en- 
ergy of the electrons, kinetic energy of the nuclei, and the 
full electron-nuclear Coulomb energy. The solute-solvent 
Hamiltonian for ET in dipolar solvents is commonly given 
by the coupling of the operator of the solute electric field 
Eq to the dipolar polarization of the solvent P 



H 0s = -E * P. 



(4) 



The bath Hamiltonian represents Gaussian statistics of 
the collective mode P with the linear response function 
X(r,r') 



*P. 



(5) 



The asterisk between the bold capital letters denotes ten- 
sor contraction (scalar product for vectors) and space in- 
tegration over the volume occupied by the solvent 



E*P 



E Pdr. 



(G) 



Assuming that the intramolecular vibrations are de- 
coupled from solvent nuclear modes allows one to cast 
the FCWD as a convolution of the vibrational, G„(w), 
and solvent, G s (w), FC densities [55|: 

/oo 
<kj'G v (uj')G 8 (u-bj' - AG/fi), (7) 
-oo 

where the diabatic equilibrium free energy gap is the sum 
of the gas-phase component AG g and difference in solva- 
tion energies AG S 



AG = AG 5 + AG S 



(8) 



In the absence of vibrational frequency change, the for- 
mer component is equal to the 0-0 transition energy in 
the gas phase. 

The FC density for each nuclear mode is given in terms 
of a broadening function g n (t) |56l l57j 

dt 



2vr 



exp[i(\ n /h-u;')t-g n (t)], (9) 



where 

9n(t) =- / -t(1 - cos zt)x"( z ) coth 

7T Jo z 



hz 

2k B T 



dz 



(10) 



(zt -sin zt)x'n(z)- 



k Jo * 

In Eq. © , A„ is the nuclear reorganization energy 

h r°° dz 
I" Jo z 



A, 



(ii) 



and x'n( z ) ls the imaginary part of the frequency- 
dependent linear response function (spectral density) cor- 
responding to the nuclear mode n (in general, many such 
modes contribute to the solvent (s) and vibrational (v) 
FC densities). 

For a set of vibrational normal modes with frequencies 
ujo and reorganization energies A g , the spectral density is 



m 



<!v(. Z ) = K^Sq^q [8(Z-U) q ) - 6(z + U) q )] 



(12) 



where S q = X q /huj q is the Huang-Rhys factor. When 
all nuclear modes are classical, g n (t) = k^TXnt 2 /h 2 and 
one reaches the classical, high temperature limit of the 
Marcus theory 



FCWDi(cj) = [4tt(A s + X v )k B T] 1/2 

(AG + A s + A.„-^w) 2 "i 



exp 



4fc B T(A s + A„) 



(13) 
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where A„ = Y] X q is the total vibrational reorganization 
energy. 

When the solvent mode is classical and the vibra- 
tions are quantized, one can use the small t expansion 
in Eq. @, valid in the limit when u q is much smaller 
than the vertical energy gap [\X V — u>'\ in Eq. ©]. With 
Uuj q /k B T ^> 1, one gets 



where 



g v (t) ~ (t 2 /2h)u v X v , 



(14) 



(15) 



is the effective vibrational frequency. With the vibra- 
tional broadenin g fu nction in the form of Eq. I|14l) the 
FCWD becomes [El |U ED 



(16) 



FCWDi(w) = [n{4:\ s k B T + 2Huj v \ v )} 1/2 

(AG + X s + X v - huj) 2 



exp 



4fcBTA s + 2Tilj v \ v 



The above equation, present in some early papers on ET 
[59l loci . is not very accurate as was pointed out by Jort- 
ner [62]. The set of equations given below, which can 
be found in work by Lax [g^, Davydov 01, and Kubo 
and Toyozawa provides a better description of the 
vibronic envelope. 

When the normal mode vibrations are represented by 
a single effective vibration with frequency defined by Eq. 
(|15f) . the vibrational FCWD is a weighted sum of reso- 
nant vibrational transitions 



G v (uj) = 2J A m 8{u-mL0 v ), 



where 



A _ _-Scothx„+rox„ T 



( S 



\sinh Xv 



(17) 



(18) 



S = \ v /fiuj Vl Xv = fu^v/2k B T, and I m (x) is the modified 
Bessel function of order to. 

The FCWD for the classical nuclear modes of the sol- 
vent is given by the expression 



Gs(u- AG/h) = h(6(AE(P n )-hu)), (19) 



where 



AE(P n ) = AG + X s - AE * SP r 



(20) 



and (. . • )i denotes an ensemble average over the fluctu- 
ations of the nuclear solvent polarization P„ coupled to 
the difference in initial and final state electric fields of 
the donor-acceptor complex, AEo = E02 — Eqi- In Eq. 
(121 H , A s stands for the solvent reorganization energy (see 
below), and <5P„ is the fluctuation of the nuclear polar- 
ization with respect to its equilibrium value. With the 



Gaussian Hamiltonian for polarization fluctuations [Eq. 
|JSJ], G s (u! — AG/h) is a Gaussian function leading to a 
total FCWD in the form of a weighted sum of Gaussians 



FCWDi(w) = [4TrX s k B T} 



-1/2 



m— — oo 



exp 



(AG + X s + mhuj v — Hu>Y 
4A s fc B T 



(21) 



When the energy of vibrational excitations is much 
greater than k B T [x v > 1 in Eq. lfT3)l] the FC envelope 
turns into a sum of Gaussians with weights given by the 
Poisson distribution 551 



■ S r ' 



.1 ' 



TO > 0. 



(22) 



III. SOLVATION THERMODYNAMICS 

Inserting a solute into a molecular solvent results in 
solvent perturbation that can roughly be split into two 
components with drastically different length scales. The 
first component is due to repulsion of the solvent from 
the solute core caused by short-range, but strong re- 
pulsive forces. This perturbation creates a local den- 
sity profile in the solvent around the solute which may 
or may not induce a polarization field acting on the so- 
lute charges. The electric field of solute charges creates 
yet another perturbation. The solute electric field is suf- 
ficiently long-ranged to induce the dipolar polarization 
P(r) in a quasi-macroscopic region of the solvent around 
the solute. Gradients of the solute field couple to the 
higher-order (quadrupolar, etc.) polarization, but this 
interaction is more short-ranged |38l l65l I66L |67j • 

The dipolar polarization is caused by alignment of the 
permanent and induced solvent dipoles along the solute 
field. This alignment occurs on two quite different time 
scales: ~ 10~ 15 s for induced dipoles and ~ 10~ n — 10~ 12 
s for permanent dipoles. Accordingly, the polarization 
field splits into a fast relaxing electronic polarization (in- 
duced dipoles, P e ) and a much slower nuclear polariza- 
tion (permanent dipoles, P„) [68|]. The electronic sol- 
vent polarization is always in equilibrium with the chang- 
ing distribution of the electronic density in the donor- 
acceptor complex. The energy conservation condition of 
the Golden Rule formula is thus imposed on the energies 
with equilibrated electronic polarization. Therefore, be- 
fore being used in the Golden Rule expression, the Hamil- 
tonian matrix should be averaged over the fast electronic 
component of the dipolar polarization [fjjj . For the 
energy Ei depending on the instantaneous configuration 
of the nuclear subsystem one gets 



-PEi 



Tr, 



-Hi/k B T 



(23) 



where Tr e i denotes the statistical average over the elec- 
tronic degrees of freedom of the solvent. Before going 
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into the details of separate calculations for electronic and 
nuclear components of the polarization, we outline the 
general formalism of polarization response to an external 
electric field. 



{') □ p M (b) 




FIG. 2: Two approaches to the calculation of the response 
function: as polarization response to an external electric field 
perturbation (a) and as correlation of polarization fluctua- 
tions near the solute hard core from which the polarization 
field is excluded (b). 



equivalent to the ensemble average (. . . )o in the presence 
of a fictitious solute which has the geometry of the real 
solute (and therefore the complete repulsion potential) 
but no partial charges. This notion provides a conve- 
nient route to the calculations of the response function 
for complex solutes. Instead of calculating the polar- 
ization in response to a non-trivial field Eo(r), one can 
calculate the correlation of polarization fluctuations in 
the presence of a fictitious solute with only the hard re- 
pulsive core of the real solute retained. The correlation 
function is then calculated with the requirement of zero 
polarization within the solute (Fig. |2Jd) 

X (r, r') = (fc B 7 1 )- 1 ( ( 5P(r) ( 5P(r'))o. (26) 

This is the essence of the approach adopted in the present 
formalism, making the response function solely deter- 
mined by the molecular structure inherent to the pure 
solvent and the short-range perturbation produced by 
the repulsive core of the solute. 



Formalism 



In the linear response approximation (LRA), the sol- 
vent polarization P(r) is a linear functional of the per- 
turbing electric field Eo (vacuum electric field of the so- 
lute for solvation): 



P(r) = X *E = f X (r,r')-E (r')dr'. 



(24) 



Here, x(r,r') is a two-rank tensor describing the non- 
local linear response of the solvent to the solute electric 
field, dot denotes tensor contraction over the common 
Cartesian projections. This function is different from di- 
electric susceptibility appearing in Maxwell's equations 
in two respects. First, x(r,r') describes the polarization 
response to the field of external charges and not to the 
total electric field E = Eo + Ep combining the exter- 
nal field with the electric field Ep created by the solvent 
polarization (x corresponds to x° of Madden and Kivel- 
son [Ifj). Second, x( r i r ') is affected by the presence 
of the solute and thus x( r i> r ') is generally not equal to 
x( r 2, r") even if ri — r' = T2 — r" . 

Equation (|24J) determines the response function in 
terms of an external electrostatic perturbation and polar- 
ization induced by it (Fig.|2t)- An alternative view of the 
response function is through the fluctuation-dissipation 
theorem which relates the response to the correlation 
function of polarization fluctuations in the solute vicinity 



X (r,r') = (fcBT)- 1 ((5P(r)5P(r')). 



(25) 



An important result of the LRA is that this correlation 
function does not depend on the long-range electrostatic 
field of the solute. The ensemble average (...) in the 
presence of the real solute with its charge distribution is 




FIG. 3: Ai for a neutral diatomic D-A (circles) and A2 for a po- 
lar diatomic D + -A~ (diamonds) for a donor- accept or complex 
represented by two contact spheres with radii Ro/a = 0.9. 
Solvent is a fluid of dipolar hard spheres of diameter a and 
dipole moment m. The simulations, reported in Ref. l36l were 
carried out at different m and at constant density pa 3 = 0.8. 
The change in solvent polarity is reflected by the dipolar den- 
sity y = (47r/9)m 2 p/fc B T. 

The applicability of the LRA to solvation of large 
donor-acceptor complexes common for ET research in 
molecular solvents is well su ppo rted by exis ting evidence 
from computer simulations [Hi E3, IZ3, E3, EJ ■ The di- 
rect consequence of the LRA are the following relations 
for the moments of the solute-solvent interaction poten- 
tial vq s : 



-k B T(v 0s ) = ((Sv 0s ) 2 ) = ((5v 0s ) 2 ) . 



(27) 



When the solute-solvent interaction is limited to the cou- 
pling of the solute charges to the solvent dipolar polar- 
ization, vq s = (*|flos|*) in Eq. @. 

The independence of the response function with re- 
spect to the solute charge is propagated into equality 
of the variance of vq s in equilibrium with fully charged 
solute, (■■■), and in equilibrium with uncharged solute, 
(. . . )n. Figure shows the results of simulations from 
Ref. |36| for a model diatomic donor-acceptor complex D- 
A in a dense solvent of hard sphere point dipoles. The 
system is designed to mimic the charge separation, D-A 



6 



— > D + -A~, and charge recombination, D + -A~ — ► D-A, 
reactions. The reorganization energies for charge sepa- 
ration, Ai = ((SvQ S ) 2 )r ) /2kj i T, and for charge recombina- 
tion, A2 = ((Svos) 2 } /2k^T , turn out to be very similar 
over a broad range of solvent polarities monitored by the 
dipolar density parameter y = (An /9)rn 2 p/k^T; p is the 
solvent number density, m is the solvent molecule per- 
manent dipole moment. 

The inhomogeneous character of the response functions 
is retained after transformation to k-space. The function 
x(k, k') then depends on two wavevectors in contrast to 
the dependence on a single wave- vector for the homoge- 
neous dielectric response. The calculation of x(k, k') is 
still a major challenge for microscopic theories of polar 
solvation. Despite some very active research in this area 
for the last 80 years since the formulation of the Born 
model for solvation of spherical ions [7r| , no microscopic 
solution applicable to solutes of arbitrary shape has been 
presented so far. A promising strategy, adopted already 
in the Born [s| and Onsager models, is to calculate 
the response functions in terms of properties of the pure 
solvent. This connection can be achieved by considering 
the polarization correlation function in the presence of 
the repulsive core of the solute [Eq. lj27fl) ]. 

The exclusion of the polarization field from the solute 
volume is provided by the Li-Kardar-Chandler approach 
[2H ITsf. in which the trajectories defining the response 
function in its path integral representation are restricted 
from entering the solute. The result of this procedure 
is an integral equation relating x(k, k') to the non-local 
susceptibility of the pure solvent X s (k) (with a single 
k-vector for the homogeneous response) and the shape 
of the solute. The equation for the response function 
is then equivalent to the Ornstein-Zernike equation for 
the solute-solvent correlation function with the Percus- 
Yevick closure for the solute-solvent direct correlation 
function ps[ . 

No general solution for x(k, k') in the Li-Kardar- 
Chandler integral equation has been obtained so far. 
One can, however, employ analytical properties of the 
response functions to obtain the solvation chemical po- 
tential [3{| 

1 f dlcrlW ~ 

-mo s = 2 J ^r E °( k ) ' *( k > k ') ' E °(- k ')- ( 28 ) 

The closed-form result for /j,q s exists when the Fourier 
transform of the electric field Eo(k) is known in analyt- 
ical functional form. This is not the case for many real 
problems, when the distribution of molecular charge is 
given from force fields or quantum calculations and the 
Fourier transform of the field is calculated numerically. 
Unfortunately, the analytical solution is given by the dif- 
ference of two large numbers almost canceling each other. 
It therefore becomes not very practical in strongly polar 
solvents because of accumulation of numerical errors. To 
facilitate numerical applications, a mean-field solution for 
x(k, k') was offered in Ref. y(| This solution eliminates 
the inhomogeneous character of the response function by 



a non-local renormalization of its transverse component: 



X(k, k') = (2^) 3 <5(k - k') [ X L (k)J L + X T (k)J T ] , (29) 

where J L = kk and J T = 1 — kk are, respectively, the 
longitudinal and transverse projections of a 2-rank tensor 
with the axial symmetry established by the direction of 
the wavevector, k = k/fc. The 6D integral of Eq. (|28|) is 
then reduced to the computationally tractable 3D inte- 
gral. 

The transverse, x T (k), and longitudinal, % L (k), pro- 
jections in Eq. I)29|l are related to corresponding compo- 
nents of the susceptibility of the pure polar solvent 



Xtr 



F .J L -E (k) 
F -J T -E (k) 



and 



In Eq. m>- Xtr 



x L (k) = x ^(fc). 

(l/3)Tr[ Xs (0)] and 
_ 2[^(0)- X f(0)] 



(30) 



(31) 



(32) 



Further, Eo(k) denotes the Fourier transform of the elec- 
tric field of the solute calculated on the volume of the 
solvent Q obtained by excluding the hard repulsive core 
of the solute from the solvent 



E (k) 



E (r)e lkr dr 



(33) 



The mean-field approximation adopted in deriving 
Eqs. (|29|I - H33|I consists of replacing a generally non- 
uniform field of the solvent within the solute by its spatial 
average Fo [Eq. i|30|) ]. The neglect of the gradients of the 
field induced by the solvent within the solute amounts to 
taking the dipolar projection of the solute field according 
to the following relation: 



where 



F = / E (r)-D r ^, 



D r = 3ff 1. 



(34) 



(35) 



is the dipolar tensor. The electric field F is a gener- 
alization of the Onsager reaction field for the case of 
non-spherical solutes with non-dipolar charge distribu- 
tion. F reduces to the Onsager field for spherical dipolar 
solutes. 

The mean-field renormalization of the transverse com- 
ponent of the response function in Eq. I|30[> resolves the 
fundamental difficulty of microscopic solvation theories 
arising from the fact that the short-range repulsive per- 
turbation caused by the solute produces a major change 
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in the polarization response functions compared to those 
of the pure solvent. For instance, a direct replacement 
of x(k) with Xs(k) in the homogeneous approximation 
(see Ref. [z^ for discussion) results in divergent behavior 
of \ s with increasing solvent dipole moment |80( . The 
divergence arises from the transverse component of the 
response ("transverse catastrophe") which has to be in- 
cluded once the dielectric cavity does not coincide with 
an equipotential surface of the solute charge distribution 
|8lj . In continuum calculations, the divergent behavior 
is eliminated by imposing boundary conditions at the di- 
electric cavity on the solution of the Poisson equation. 

Although the problem with the transverse response 
has long been recognized in the literature 0, 182^ . 
many microscopic formulations of solvation thermody- 
namics and dynamics have avoided the problem by ne- 
glecting the transverse response [2^, [83l |§4| which is also 
neglected in some continuum calculations, e.g. the Gen- 
eralized Born approximation j8^. Equations (|29[) -H33 [I 
provide a general solution of the problem which ag rees 
well with available simulations of polar solvation [35[ |3|J 
and experiment on solvation dynamics |39j . The formal- 
ism is based on the homogeneous solvent susceptibility 
as input and, once the susceptibility is defined from com- 
puter experiment or liquid-state theories, can be applied 
to solvation in an arbitrary isotropic dielectric. 




FIG. 4: Components of the induced dipole moment at sol- 
vent molecule j. p° is produced by the external electric field 
Eo(rj); pf is produced by the reaction field induced in the 
solvent by the permanent dipole nij. 



B. Polarization structure factors 

The dipole moment at a given molecule j in a polar- 
polarizable solvent is a sum of the permanent dipole nij 
and the induced dipole pj 

Vj = m j+Pj- ( 36 ) 

The total induced dipole then splits into p° created by 
the external electric field Eo(rj) and pf induced by the 



reaction field (superscript "R" ) caused by the dipole rn, 
itself (Fig. HJ: 

Pj = P ■ + Pf (37) 

The reaction field caused by the dipole nx, relaxes on the 
time-scale of translational-rotational motion of molecule 
j. Therefore, the induced dipole pf , which follows adi- 
abatically the reaction field, should be attributed j8|| 
to the slow nuclear polarization of the solvent P„. In 
contrast, the component p°, following adiabatically the 
external field, is attributed to the the fast solvent polar- 
ization P e . The sum of the permanent dipole ixij and the 
induced dipole pf makes the effective condensed-phase 
dipole H3 

m'j = m, + pf = m'ej , (38) 

where e., is the unit vector along the direction of rxij. 
The dipole moment m! in principle depends on the in- 
stantaneous configuration of the liquid. However, we will 
not consider fluctuations of m' here an(L following self- 
consistent models of polarizable liquids [87|, will replace 
m' with its statistical average value. 

The attribution of the electronic polarization in equi- 
librium with the electric field of the permanent dipolcs 
to the nuclear (slow) polarization of the solvent is an 
essential part of the Pekar partitioning of the solvent po- 
larization into fast and slow components l88l |89l . Other 
partitioning schemes have been proposed 90], but they 
all lead to the same value of the solvation energy when 
correctly implemented |9lT |. Computer simulation proto- 
cols in which the induced polarization is self-consistently 
adjusted to the instantaneous nuclear configuration pro- 
vide direct access to the slow polarization in Pekar's defi- 
nition |92(. Self-consistent simulations of polarizable sol- 
vents are used here to test the analytical procedure em- 
ployed for the response functions of the nuclear polariza- 
tion (Sec. HQ} . 

The total dipolar response function of the homoge- 
neous solvent is a 2-rank tensor describing correlations 
of dipole moments fij: 

X.(k) = (fi/Sl) (j^ MiMAe*"" ^ , (39) 

where = — 17. and brackets refer to an ensemble av- 
erage. Because of the isotropic symmetry of the solvent, 
g|k) splits into longitudinal and transverse components 

Xs (k)=Xs(k)J L +X^(k)J T . (40) 

It is convenient to factor the response function into the 
effective density of dipoles y c s, which is mostly affected 
by the magnitude of the solvent dipole, and the struc- 
ture factor, which reflects dipolar correlations and can 
be expressed through angular projections of the pair dis- 
tribution function 36] 

Xs W = ^[S L (k)J L + ^(k)J T }. (41) 
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The structure factors S L ' T (k) (Fig. |SJ are denned based 
on the unit vectors \ij = fjbj/ fij in the direction of the 
respective total dipole moments 



dielectric constant e, 



(42) 



The effective dipole density in Eq. (|4i() is 

2/eff = V P + (47r/3)pa, y p = (4ir/9)p(m') 2 /k B T, (43) 

where a is the dipolar polarizability. Only the perma- 
nent dipole moment is renormalized by the mean field of 
the solvent in the above e quat ion, which corresponds to 
Wertheim's 1-RPT theory [93| (2-RPT theory renormal- 
izes the polarizability a to a', but the 1-RPT version of 
the theory is in better agreement with simulations |l9j]h 
The nuclear response function reflects correlated ori- 
entations and positions of dipoles m^ : 



X„(k) = (/3/^)(^m>^ k - r 



(44) 



Similarly to Eq. 1|41J) . Xn(k) can be separated into the 
longitudinal and transverse components 



Xn (V = 3 -^[sf;(k)J L + s? l (k)j T ] 



(45) 



The nuclear structure factors are defined by Eq. H42|) . in 
which the unit vectors (ij are replaced by the unit vectors 
Bj [Eq. OEM- 




FIG. 5: Longitudinal (L) and transverse (T) polarization 
structure factors calculated by using the PPSF with the pa- 
rameters of ambient water. The solid lines refer to the total 
structure factors S L ' T (k), the dashed lines refer to the nuclear 
structure factors S^' T (k). 

The k = values of the structure factors are related to 
the macroscopic dielectric properties of the solvent. The 
total polarization response is defined through the static 



S L {0) 

s T (o) 



e s -l 
3e s y e ff ' 
e s -l 
3y c s 



(46) 



The nuclear structure factors depend, in addition, on the 
high-frequency dielectric constant Coo 1123 



3%,' 
3y P 



where 



Co = 1/coo - l/e s 



(47) 



(48) 



is the Pekar factor. 

Both S^ T (k) and S L - T (k) tend to unity at k — > oo. 
This limit is the result of the point multipole approxi- 
mation for the intramolecular charge distribution within 
the solvent molecules. In contrast, charge-charge struc- 
ture factors defined on interaction-site models of liquids 
decay to zero at k -> oo (Refs. |H HI H3). Thc re- 
gion of k- values where this distinction becomes important 
is, however, insignificant for the calculation of solvation 
thermodynamics (see below). The nuclear and the total 
structure factors differ in the range of small k- values and 
around the longitudinal peak as a result of the influence 
of the high-frequency dielectric constant of the solvent 
(Fig. EJ. The effect of on the longitudinal peak is 
insignificant for the calculation of the reorganization en- 
ergy. Therefore, it is the range of small fc-values and, in 
addition, the dependence of the liquid-state dipole mo- 
ment to' on the solvent polarizability, that ultimately 
determine the variation of the solvent reorganization en- 
ergy with the solvent high-frequency dielectric constant 
Coo (see below). 



C. ET thermodynamics 

The solvation thermodynamics of ET is determined by 
the solvent reorganization energy and the solvent compo- 
nent of the free energy gap. They are defined in terms of 
the nuclear and total response functions by the following 
relations 

Xs= 2j T2^ AE ° (k) ' X " (k ' k ' } ' AE °(- k ') ( 49 ) 



and 



dkdk' 



AE (k)-x(k,k')-E (-k'). (50) 



In Eqs. EH and AE (k) = E 02 (k) - E i(k) and 

Eo(k) = (E 02 (k) +E i(k))/2; E 0l (k) are the Fourier 
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transforms of the solute electric field in the initial (i = 1) 
and final (i = 2) ET states taken over the volume O 
occupied by the solvent [Eq. (J33}]. 

The mean-field solution for the response functions [Eq. 
(|29[) ] splits both the solvent reorganization energy and 
the free energy gap into their corresponding longitudinal 
and transverse components: 



and 



AG, = AGf + AG^ 



(51) 



(52) 



Each projection is obtained as a k-integral with the cor- 
responding polarization structure factor. For the "T" 
projections one gets 



A 



3y p #(0) 



dk 



and 



8tt g Kn J (2tt) 3 
ag t = 3y cS S L (0) 



A£ T (k) Sl{k) (53) 



rfk 



8tt g K J (2^) 3 
l^o T 2 (k)| 2 -|^ T i(k)| 2 ] S T {k). 
In Eqs. and 

S jf„ = (l/3) [#(0) + 2S£(Q)] 

and 

to = (1/3) [S L (0) + 2S T (0)] 



(54) 



(55) 



(56) 



are the nuclear and total Kirkwood factors, respectively. 

The longitudinal components of free energies, Af and 
AGg , include both the longitudinal and transverse pro- 
jections of the solute field: 



3y P 



dk 



8tt J (2?r) 



(57) 



and 
AG^ 



dk 



8tt J (2tt) 3 
In Eqs. and 0, 



^(k)-^))^). (58) 



fif (k) = |A£#(k) 



and 



fnlAE ° (k)l AF .J-.AE„(k) 

(59) 



£f(k) = | J B L i (k)| 2 -/ s | J B T l (k)| 



> Foi • J L • Eoi(k) 
F 0j -J T -E 0l (k)' 



(60) 



The longitudinal and transverse components of the elec- 
trostatic energy density in Eqs. are defined as 



A£ L ' T (k)| 2 = AEo(k) • 3 L ' T ■ AEo(-k), 



E^ T (k)\* = E 0l (k) • J^' • E w (-k) 



L,T 



(61) 



The effective fields £^(k) and £f s (k) depend on the sym- 
metry of the charge distribution within the solute analo- 
gously to the result of imposing the boundary conditions 
on the solution of the Poisson equation in continuum elec- 
trostatics. 

The electric field Foi in Eq. (|60|l is a generalization of 
the Onsager reaction cavity field [73 to the case of so- 
lutes of non-spherical shape and non-point-dipole charge 
distribution. This field is obtained by summing up a con- 
tinuous distribution of dipolar electric fields induced by 
the solute in the solvent volume: 



f E w (r).D r *, 



(62) 



where D r is given by Eq. (J3SJ). Also, AF in Eq. (JSSJl is 
AF = F 2 — Foi. F i becomes the standard Onsager 
reaction field for a point dipole at the center of a spherical 
cavity. Finally, in Eqs. JS^I and (|dT7|| . 



2(e s -l) 



2e, 



1 



(63) 



is the usual Onsager polarity parameter [23] and the cor- 
responding polarity parameter for the nuclear polariza- 
tion is 



2(eoo£s - 1) 



1 



(64) 



IV. CALCULATION PROCEDURE 

The formalism outlined above is realized in a compu- 
tational algorithm sketched in Figure El It includes two 
branches, one is for the solvent part of the calculation and 
another is for the solute part. The two parts are com- 
bined together in the integration over the inverted space, 
which yields the reorganization energy (A s ) and the to- 
tal free energy of nuclear plus electronic solvation (AG S ). 
We start with describing the solvent branch followed by 
the outline of the solute part. 



A. Solvent 

The calculation of the structure factors in the sol- 
vent branch in Fig. requires a set of experimental in- 
put parameters: m (gas-phase dipole moment), a (gas- 
phase dipolar polarizability) , (high- frequency dielec- 
tric constant), e s (static dielectric constant), and a (ef- 
fective hard sphere diameter of the solvent molecules). 
The hard sphere diameter is obtained from the exper- 
imental compressibility of the solvent by fitting it to 
the compressibility found from t he g eneralized van der 
Waals (vdW) equation of state 95]. Based on these 
parameters, an analytical procedure has been recently 
proposed to calculate S L ' T (k) (3(|. This parameteriza- 
tion, called parametrized polarization structure factors 
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Solvent parameters: 
m, a, p, 6 , a 



Solute parameters: 

X 0k' ^Ok' Z 0k' %k' a 0k 



S L (k), S T (k) E(k) 



k-integration 



A., AG 



FIG. 6: Diagram of the calculation algorithm. Solvent pa- 
rameters include: m (gas-phase dipole moment), a (gas- 
phase dipolar polarizability) , too (high-frequency dielectric 
constant), e s (static dielectric constant), and a (effective hard 
sphere diameter of the solvent molecules). Parameters xok, 
yok, zok stand for Cartesian coordinated of the solute atoms, 
qok are partial charges, and anfc are atomic vdW radii. 



(PPSF), makes use of the analytical solution of the mean- 
spherical approximation (MSA) for dipolar fluids [96| . 
The MSA solution gives S L,T (k) in terms of the Baxter 
function Q(ka, rj) appearing as solution of Percus-Yevick 
integral equations for hard sphere fluids [97j 



S(ka,rj) = \Q(k<r,T))\~ 



where 



■ikat 



g(fccr,r;) =1 - 1277 / e" 
Jo 

[a( V )(t 2 -l)/2-b( V )(t-l)]dt 



(65) 



(66) 



and a{ri) = (1 + 2? ? )/(l - 7;) 2 , 6(r?) = -3?j/2(l - t?) 2 . 
For a fluid of hard sphere molecules, r\ = (ir/6)pcr 3 is 
the packing density, equal to the ratio of the volume of 
the solvent molecules to the volume of the liquid. In the 
MSA, the S L ' T (k) are obtained by setting 77 — 2£ for 
S L (k) and 7/ = -£ for S T {k) in Eq. @5J- Here, f is the 
MSA polarity parameter which can be related either to 
th e dip olar density y c g or to the static dielectric constant 

Two problems arise when dealing with the reorganiza- 
tion energy calculations using the polarization structure 
factors from the MSA. First, one needs a general proce- 
dure which would provide the nuclear structure factors 
Sn' T (k) in polarizable solvents in contrast to total struc- 
ture factors S L ' T (k) given by the MSA solution. Such a 
formalism should thus exclude (quantum) fluctuations of 
the induced solvent dipoles p9 which are not included in 
the nuclear polarization field (Fig. 0J. Second, the MSA 
does not give a consistent description of the dielectric 
properties of polar solvents, i.e. the polarity parameters £ 
calculated from y e g and e s are quite different. The PPSF 
procedure goes around the second problem by consider- 
ing y c g and e s as two independent input parameters used 
to calculate S L ' T (k). A convenient way to introduce the 
two-parameter scheme is to specify two separate polarity 



parameters which are obtained from the longitudinal and 
transverse structure factors at k = 0: 



(l-2£ 



L\4 



(l + 4^) 2 
(1-2H 2 



=s L (o), 

=6^(0). 



(67) 




FIG. 7: Nuclear longitudinal (L) and transverse (T) struc- 
ture factors from the PPSF (dashed lines) and MC simu- 
lations (solid lines). MC simulations are carried out for a 
fluid of 1372 hard spheres with permanent dipole m, diame- 
ter a, polarizability a, and density p: (m*) 2 — /3m 2 /a 3 — 1.0, 
a* — a/a :i = 0.06, pa 3 = 0.8. The dielectric properties from 
the simulations are: e s — 21.4, y e fi = 1.57, and y p — 1.54; 
too = 1.75 is obtained from the Clausius-Mossotti equation. 

Separate definitions of £ and £ T in terms of S L (0) 
and 5 T (0) [Eq. gjjj] allows us to incorporate contribu- 
tions to macroscopic dielectric properties which are not 
present in the model of dipolar HS fluids. Specifically, 
the magnitude of parameter y e ft, calculated according 
to Wertheim's 1-RPT algorithm 93], defines the solvent 
dipolar strength which strongly affects the dielectric con- 
stant. However, e s also depends on such factors as sol- 
vent quadrupolar moment |87j|. solvent non-sphericity, 
etc. The influence of these factors is incorporated into 
S L ' T (0) through the dielectric constant. Similarly, the 
polarity parameters and £j are calculated from Eq. 
(E3 with S L ' T (0) replaced by S£' T (0) taken from Eq. 

Dipolar projections of the structure factors of molecu- 
lar liquids modeled by site-site inter actio n potentials have 
been studied previously HI HI GjIJ . The PPSF pro- 
cedure has also been tested against MC simulations of 
dipolar hard sphere fluids [3(j. However, the structure 
factors arising from the nuclear polarization as well as 
the applicability of the PPSF to non-spherical molecules 
with site-site potentials have not been previously tested. 
This is the aim of the Monte Carlo (MC) and MD sim- 
ulations carried out in this study. The details of the 
simulation protocol are given in Appendix ^ and here 
we focus only on the results. 

Figure \7\ shows the comparison of the transverse and 
longitudinal components of the nuclear structure factors 
calculated from the PPSF and from MC simulations. 
The MC simulations (dashed lines in Fig. [JJ have been 
performed on a fluid of 1372 polarizable dipolar hard 
spheres characterized by dipole moment m, diameter a, 
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and isotropic polarizability a ((to*) 2 = (3m 2 /er 3 = 1.0, 
a* = a/a 3 = 0.06, Appendix |A"|) . Since the simulation 
protocol generates the induced polarization in equilib- 
rium with the nuclear configuration of the solvent [lfij . 
the generated ensamble yields the nuclear polarization in 
the Pekar partitioning 88]. 

The PPSF nuclear structure factors are calculated by 
the relations: 

S%{k) = \Q(ka, -&)]-* (68) 

and 

S%(k) = \Q{ K ka,2ti)\- 2 . (69) 

In Eq. (|69|l . K — 0.95 is an empirical parameter intro- 
duced for a better agreement between the PPSF and MC 
simulations of non-polarizable dipolar fluids j93|. The 
simulations and the PPSF agree well in the entire range 
of solvent polarizabilities a* = a /a 3 = 0.01 — 0.08 stud- 
ied by simulations 92]. 




k.A' 



FIG. 8: Upper panel: longitudinal (1), k 2 {£f(k))} & , and 
transverse (2), k 2 ((AEq (k)) 2 )^ components of the electro- 
static energy density of complex 1 entering the fc-integral in 
Eqs. 1571 and 15311 . respectively. (...)£ denotes the average 
over the orientations of the wavevector k. Lower panel: lon- 
gitudinal (L) and transverse (T) structure factors for TIP3P 
water at 298 K. The solid lines refer to the results of MD 
simulations. Dashed lines indicate the results of PPSF calcu- 
lations with the input parameters corresponding to the TIP3P 
force field (Table [VJ m = 2.35 D, e s = 95.4, e x = 1.0) and 
a — 2.87 A. The dash-dotted lines refer to the nuclear struc- 
ture factors of ambient water. The graphs in the upper and 
lower panels are plotted against the same scale of k- values to 
indicate that the details of the structure factors beyond ap- 
proximately ka ~ 7r are insignificant for the calculation of the 
reorganization energy. 

The MSA solution in Eq. i|65|) was derived for a model 
liquid of dipolar hard spheres. The parameterization in- 
troduced by the PPSF suggests to use the experimental 



e s to accommodate empirically the features which are not 
included in the MSA solution. Two factors, often present 
in real polar solvents, molecular quadrupoles and non- 
sphericity, are expected to affect significantly the form of 
the structure factors. Therefore, we have performed MD 
simulations for t wo solvents with wel l-developed force 
fields, water |101| and acetonitrile |102| . Water is a rela- 
tively sym metric molecule with a very large quadrupole 
moment Q [l03| ((Q*) 2 = PQ 2 /cr 5 = 1.1) among com- 
monly used molecular solvents. On the other hand, ace- 
tonitrile has a small quadrupole moment ((Q*) 2 — 0.13), 
but the molecule is very non-spherical with the aspect ra- 
tio ~ 3. Therefore, these two extreme cases may provide 
a good test of the ability of the PPSF to incorporate the 
complications related to molecular specifics of the sol- 
vents in terms of their macroscopic dielectric constants. 

Figure [S] (lower panel) shows the comparison of the 
simulation results for TIP3P water to the PPSF. A 
slightly wrong positioning of the longitudinal peak may 
be related to a different hard sphere diameter of TIP3P 
water (see TablelVlin Appendix lA"|l compared to the hard 
sphere diameter of water at ambient conditions used in 
scaling wavevectors in Figure [S] A downward scaling of 
a by just 5% results in a very good match between calcu- 
lated and simulated structure factors. As expected, the 
steric effects of packing the solvent molecules in dense 
liquids is the main factor determining the position of the 
longitudinal peak. This indeed is seen in Fig. for simu- 
lations of acetonitrile. The effective hard sphere diameter 
obtained from solvent compressibility does not accommo- 
date the fact that linear dipoles tend to pack side-to-side 
pointing in opposite directions. The longitudinal thus 
peak effectively reflects a lower molecular diameter. The 
preferential opposite orientation of the dipoles leads to 
a low Kirkwood factor and the dielectric constant much 
lower than one would expect for a dipolar solvent with 
such large dipole moment (4.12 D f or th e force field by 
Edwards, Madden, and McDonald 102]). As a result, 
the transverse structure factor does not change with k as 
much as it does for hard sphere dipolar liquids (cf. Figs. 
□ and OH to Fig. 0. As is seen, the PPSF with e s from 
MD simulations accommodates this feature of the solvent 
quite well. 

Figure |S| compares on the common scale the k- 
dependence of the longitudinal and transverse compo- 
nents of the electrostatic energy density of complex 1, 
k 2 (£f(k))^ and fc 2 (|A^(k)| 2 ) £ , with the longitudinal 
and transverse components of the polarization structure 
factors ((■■■)£ refers to the average over the orientations 
of the wavevector k). This comparison shows that de- 
tails of the molecular structure of the polar solvent af- 
fecting the range of A:- values beyond the limit of k ~ ir/a 
are insignificant for the calculation of the reorganiza- 
tion energy and the free energy gap. Therefore, the 
discrepancies in the position of the longitudinal peak 
between the simulations and the PPSF do not notice- 
ably affect the results of calculations. This statement 
also applies to the range of k- values (k > 2ir/l s , where 
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l s is the characteristic distance between partial charges 
within the solvent molecule) at which the multipolar ap- 
proximation for the charge distribution within the sol- 
vent molecules breaks down. The charge-charge struc- 
ture factor s calculated on site -site interaction potentials 

EE HI: [ToE [IHE then deca y t0 zer0 instead 
of approaching the unity limit (S L ' T (k) — > 1 at k — > oo) 
of multipolar approximations [iM llOOl Il06j . The range 
of fc-values where the inaccuracy of the multipolar ap- 
proximation becomes significant lays beyond the range 
of small fc-values affecting the calculation of thermody- 
namic properties unless the solute is much smaller than 
the solvent. 
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FIG. 9: Longitudinal (L) and transverse (T) polarization 
structure factors of acetonitrile at 298 K. The solid lines re- 
fer to the results of MD simulations. Dashed lines indicate 
the results of PPSF calculations with the input parameters 
corresponding to ACN3 (Table El m = 4.12 D, e 3 = 29.6, 
too = 1.0) and a = 4.14 A. The dash-doted lines indicate the 
nuclear structure factors S^' T (k) from the PPSF with the 
parameters of ambient acetonitrile: m — 3.9 D, e s = 35.9, 
too = 1.8, a = 4.48 A 3 , r) = 0.424, a = 4.14 A. 



B. Solute 



the solute and Mq partial charges qok to be used in Eq. 
dZOIl (indicated as x fc, Vok, z k, qok in Fig.[BJ|. 

The infinite-space Fourier transform of the Coulomb 
electric field [Eq. is numerically divergent |36]. This 
numerical problem is obviated by splitting the region of 
integration into the inner part between the SAS and a 
cutoff sphere and the region outside the cutoff sphere. 
The Fourier transform within the sphere is calculated 
numer ically by the Fast Fourier Transform (FFT) tech- 
nique 108J on a cube with the center at the geometrical 
center of the DSA complex 

No 

r c = N^J2rok- (71) 

k=l 

The length of the cube is chosen by multiplying the max- 
imum extension of the molecule measured from r c by a 
factor of 9. This choice yields a sufficiently small incre- 
ment of the k-grid necessary for the inverted-space inte- 
gration and, at the same time, avoids numerical errors 
arising from artificial periodicity imposed by a finite-size 
numerical FFT technique. The FFT calculation was done 
on a grid of dimension 256 x 256 x 256 and the step of 0.5 
A. Calculations on complex 1 involved 143 atoms holding 
partial charges q Qk . The charge shifts (Aq k = q$ k — qQ k ) 
and coordinates used in the solvent reorganization and 
free energy calculations are the same as those reported 
in Ref. 40. The individual (i.e., initial and final state) 
charges used in the reaction free energy calculations are 
also taken from Ref. 0. The field Eni(k) obtained by 
combining the numerical and analytical parts is used to 
calculate the longitudinal and transverse components of 
the electrostatic energy density in Eqs. (|60|l and l|61|l . 
These components are then used in the /c-integrals with 
the polarization structure factors (Eqs. i|53[l ~ H58[l : also 
see Fig. EJ. 



The solute branch of the calculation algorithm (Fig. 
EJ) consists of the numerical calculation of the Fourier 
transform of the electric field outside the solute placed in 
the vacuum. The direct-space electric fields in the initial 
and final states of the solute are given by a superposition 
of electric fields produced by partial charges q l Qk 

Mo _ 

^"■Wf^ (70) 

where the sum runs over Mq partial charges localized on 
solute atoms. The field Eoj(r) is Fourier transformed in 
the region O accessible to the solvent molecules [Eq. (|33|) ] . 
The region ft is generated by assigning vdW radii to all 
atoms of the solute and then adding the hard sphere ra- 
dius a/2 of the solvent (<r = 2.87 A for water and 4.14 
A for acetonitrile). This creates the solvent-accessible 
surface (SAS). The definition of the solute field thus re- 
quires atomic coordinates and vdW radii of Nq atoms of 



V. RESULTS AND COMPARISON TO 
EXPERIMENT 

A. Solvent reorganization energy 

The solvent reorganization energy of complex 1 was 
previously obtained from MD simulations of this com- 
plex in TIP3P water ^(j- The permanent dipole moment 
in this force field is enhanced from the vacuum dipole of 
1.87 D to 2.35 D to account for water polarizability. This 
results in a dielectric constant of e s — 97.5 from our sim- 
ulations, which agrees well with e s = 97.0 found in the 
literature [l09j. Table [I] lists the results of calculations of 
the reorganization energy with structure factors from the 
PPSF (column 5) and from MD simulations (column 7). 
The density of the solvent in the NVT simulations was 
adjusted at each temperature in order to reprodu ce th e 
expansivity a p = 2.96 x 10~ 4 K" 1 of TIP3P water [TTo| . 
The temperature derivative of the reorganization energy 
thus gives the constant-pressure reorganization entropy 
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TABLE I: Reorganization energy (kcal/mol) of 1 in water. 
All calculations except those in the last column are done with 
too = 1.0. 



T, K 



A 



Ac 



"aT" 



\ J \ J 



288 0.4110 6.44 107.7 64.93 39.60 64.35 1.55 4.37 45.88 

293 0.4104 6.32 102.1 64.52 39.58 64.26 1.45 4.33 45.47 

298 0.4098 6.21 97.5 64.11 39.56 63.93 1.38 4.25 45.07 

303 0.4092 6.09 96.0 63.67 39.55 63.52 1.42 4.13 44.68 

308 0.4085 5.99 93.7 63.25 39.54 62.98 1.22 4.14 44.30 

"Packing fraction calculated with a = 2.87 A and the isobaric 
expansion coefficient a p = 2.96 X 10~ 4 K _1 . 

"TIP3P water has the permanent dipole of 2.35 D scaled up from 
the vacuum dipole moment of water, 1.83 D, to account for the 
mean-field effect of the induced dipoles. 

c Calculated from MD simulations as described in Appendix A. 

^Calculations with the PPSF structure factors with eoo = 1.0 
and e B from MD simulations. A p stands for the reorganization 
energy arising from the interaction between the solute electric field 
and the solvent dipoles, \ q comes from the interaction between the 
solute field gradient and solvent quadrupoles, X pq is the mixed term 
from correlated fluctuations of dipoles and quadrupoles on different 
solvent molecules, see Eq. 1731 . 

e Continuum limit S„' T (fc) = S L ' T (0) at = 1.0 and e s from 
MD simulations. 

■'Calculations with the structure factors from MD simulations. 

Calculations with the PPSF structure factors with the solvent 
parameters of ambient water. 



corresponding to conditions normally employed in exper- 
iment, 



s x = -(d\ s /dT)p. 



(72) 



Overall, there is an exceptionally good agreement be- 
tween the reorganization energies calculated by using the 
structure factors from PPSF and MD simulations. This 
is not surprising in view of the very good agreement be- 
tween the two sets of structure factors shown in Fig. |SJ 

The PPSF result at 298 K, A s = 64.11 kcal/mol, also 
compares well with the direct calculation of the reorga- 
nization energy from MD simulations, where the value of 
60.9 kcal/mol was reported |4fJ. The electrostatic forces 
in those simulations were cut off at distances greater than 
10.1 A. The cutoff is expected to lower the reorganization 
energy compared to that of an infinite system. In order 
to estimate the effect of the interaction cutoff, we have 
calculated the reorganization energy for a fictitious so- 
lute with the distance 10.1 A added to the radius of each 
atom exposed to the solvent. This contribution amounts 
to 7.1 kcal/mol. Column 6 in Table |I] shows the results 
of calculations when the fc-dependent polarization struc- 
ture factors are replaced by their k = values. The gap 
in A s values between columns 5 and 6 thus quantifies the 
contribution of the non-local part of solvent response to 
the reorganization energy. The last (10) column in Table 
[I] shows the PPSF calculations using parameters of am- 
bient water. In these calculations, the gas phase dipole 
moment m = 1.87 D is renormalized by the polarizability 
effect to give ml = 2 A3 D (Wertheim's 1-RPT formal- 
ism 

|Hl23)- Despite this renormalization, A s in this 



TABLE II: Reorganization energy (kcal/mol) and reorganiza- 
tion entropy (e.u., cal K _1 mol _1 ) at T=298 K of complex 1 
(experimental parameters for ambient water) calculated with 
the PPSF for the structure factors. 



Coo 


e s 


X s a 


S x a 


A s " 


5/ 


A s c 


K d 


1.0 


78.0 


52.97 


46.50 


39.49 


-8.14 


81.12 


43.64 


1.2 


78.0 


48.84 


52.17 


32.84 


-4.85 


69.61 


37.44 


1.4 


78.0 


46.41 


61.52 


28.11 


-2.88 


61.39 


33.01 


1.6 


78.0 


45.32 


70.66 


24.52 


-1.60 


55.23 


29.69 


1.8 C 


78.00 


45.15 


80.01 


21.74 


-0.71 


50.44 


27.11 


2.0 


78.0 


45.68 


90.12 


19.52 


-0.15 


46.60 


26.05 


1.0 / 97.5 9 


64.11" 84.13* 


39.56" 


2.96" 


81.44 


43.79 



"NRFT with the PPSF for ambient water with varying polariz- 
ability a. 

"Calculated with S L (k) = S L (0) and S T (k) = S T (0). 
c DelPhi calculation with the vdW cavity. 
d DelPhi calculation with the solvent-accessible cavity. 
Tarameters corresponding to water at ambient conditions. 
^TIP3P water. 

9 From present MD simulations. This value is in good agreement 
with e B = 97.0 reported in the literature |l0flj| . 

"Calculated for TIP3P water using the PPSF. 

t de s /dT = —0.654 K from MD simulations is used; this value 
turns to be higher than experimental de s /dT = —0.398 K _1 . 



calculation is substantially (~ 30 %) smaller than in the 
calculations using parameters of TIP3P water. TIP3P 
water thus appears to produce stronger solvation than 
ambient water. 

Tablc|l]also presents two components of the solvent re- 
organization energy produced by solvent quadrupoles: X q 
is the second cumulant of the coupling of the solute elec- 
tric field gradient to solvent quadrupole moment |23, E3 
whereas X pq is a mixed term arising from correlated fluc- 
tuations of dipoles an d quadru poles positioned at differ- 
ent solvent molecules 1671 111]. The resulting solvent re- 
organization energy is the sum of the dipolar component 
A p and two quadrupolar components: 



At, + A D0 + A„ 



(73) 



The problem of quadrupolar solvent re org aniz ation 
has recently attracted much attention [65l l66l l67l Ill2j 
in connection with new experimental data showing ap- 
preciable solvent reorgan ization in non-dipolar solvents 
|lll ITT1 [Till ITU ITl% However, the components X q 
and Xp q constitute only a small fraction of the overall 
reorganization energy despite a relatively high reduced 
quadrupole of water, (3Q 2 /a b = 1.1 (cf. to (Q*) 2 = 0.13 
of acetonitrile). For the rest of the paper we will therefore 
assume 



As — A,, 



(74) 



We note that the value A s = 69.7 kcal/mol calculated for 
TIP3P water with the account of water quadrupoles is in 
remarkable agreement with A s = 68 kcal/mol obtained 
by correcting the simulated values |40l | by the finite-size 
cutoff effects. 

The dependence of X s and the reorganization entropy 
on are given in Table [H] In these calculations, the 
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vacuum dipolc moment of water, 1.83 D, was held con- 
stant along with the total dielectric constant e s = 78.0. 
The change in was achieved by varying the polariz- 
ability a according to the Clausius-Mossotti equation 



1 



8rja/a 3 



(75) 



where r\ = (tt /6)pa is the solvent packing fraction. 

Two drastically different predictions for the effect of 
solvent polarizability on X s can be found in the liter- 
ature. The classical Marcus two-sphere model pj pre- 
dicts a drop of X s by about a factor of 0.6 when going 
from = 1.0 to = 1.8. On the other hand, simu- 
lations using non-polarizable and polarizable versions of 
the water force field predict almo st no dependence of A s 
on solvent polarizability 0, Ill8| . The actual situation 
is in between of the two extremes. The reorganization 
energy does drop with increasing eoo, but not as much 
as is predicted by continuum models 19]. On the other 
hand, the change is sufficient to make simulations based 
on non-polarizable solvent models unreliable. 



i — 1 — i — 1 — i — 1 — i — 1 — ; 

cont./vdW 

ambient water 




FIG. 10: A„ VS 6 co calculated for complex 1 by using the non- 
local polarization response theory (NRFT, solid lines). The 
dashed lines refer to the numerical solution of the Poisson 
equation with the vdW (cont./vdW) and SAS (cont./SAS) 
cavities. The diamond and square indicate TIP3P and ambi- 
ent water, respectively. 

The situation for the dependence of X s on is illus- 
trated in Fig. 1101 where continuum results for comp lex 1 
obtained with the DelPhi Poisson-Boltzmann solver |l!9j 
are compared to the calculations within the NRFT. The 
dielectric calculations with the vdW dielectric cavity (de- 
noted "cont./vdW" in Fig.EJ show a substantial drop of 
A s with eoo. The dependence on £ m is much weaker in the 
NRFT (see also Table ITT|) . The weak dependence of A s 
on eoo is the result of the cancellation of two competing 
factors: the decrease of the longitudinal structure factor 
in the range of small fc-values with increasing f m (Fig.[5J) 
compensated by an increase in y p due to higher solvent 
dipole m' in more polarizable solvents. We note that this 
cancellation is strongly affected by the fc-dependence of 
the polarization structure factors in the range of small 
k- values contributing to the fc-integral and cannot be re- 
duced to the cancellation of the y p factor in A s [Eqs. l(5l5|) 
and (|57|l ] with y p in the denominator in Eq. (|47|l . result- 
ing in the Pekar factor of continuum electrostatics. 



The continuum limit of the NRFT is obtained when 
the dependence on the wavevector k is neglected in the 
solvent structure factors and one assumes S L > T (k) ~ 
S L > T (0) and S%> T (k) a S%' T (0). When this assumption 
is incorporated in the microscopic calculations (marked 
NRFT/S'(0) in Fig. I10[l . the resultant reorganization en- 
ergy gains the strong dependence on characteristic of 
continuum theories. The continuum limit of the micro- 
scopic theory corresponds, however, to the dielectric cav- 
ity coinciding with the SAS. The corresponding DelPhi 
calculation (marked cont./SAS in Fig. llOfl indeed goes 
parallel with the continuum limit of the NRFT. The dis- 
tinction between these two results arises from the mean- 
field approximation used in the NRFT formulation and 
different handling of the polarizability effects in the two 
formulations (additive in the continuum and non-additive 
in the microscopic formulation |92(). Note that the mean- 
field approximation is more accurate, when compared to 
the exact solution of the Li-Kardar-Chandler equation, 
in the full microscopic formulation than in its continuum 
limit [3(j . The exact formulation of the theory, which 
does not involve the mean-field approximation, gives the 
solution of the Poisson equation in its continuum limit. 

The numerical values for the reorganization energies 
shown in Fig. ^| are given in Table Hi) The comparison 
between the microscopic and continuum calculations is 
instructive. At — 1, X s from the vdW continuum 
is much higher than the microscopic calculation, while 
A s from the SAS continuum is close to the microscopic 
result. With increasing eoo, on the other hand, X s from 
the vdW continuum falls down almost to the microscopic 
value. The continuum calculation with the vdW cavity 
may thus appear in a reasonable accord with microscopic 
calculations or experiment due to the mutual cancellation 
of errors. 

Along with reorganization energies, Table [H] lists re- 
organization entropies S\ [Eq. (|72l) ]. Note that S\ ob- 
tained from the PPSF calibrated on TIP3P water is in a 
reasonable agreement with the corresponding value ob- 
tained with the structure factors from MD simulations: 
84.1 e.u. and 69.9 e.u., respectively. The dielectric con- 
tinuum calculation gives the wrong sign for the entropy 
in accord with previous reports |12l l29j ] . Also the mag- 
nitude of S\ is substantially higher in the microscopic 
theory than in the continuum calculation (cf. columns 4 
and 6 in Table HTJ) . 

A similar trend is seen for the reaction free energy gap 
( Table ITTTf) for which the reaction entropy is defined as 



AS S = - (dAG s /dT) P . 



(76) 



Although the sign of AS S is correct in the continuum 
calculations, the entropy magnitude is much lower than 
in the NRFT, similar to a previous report for a differ- 
ent ET system |12(, where AS S was experimentally ob- 
tained from temperature dependent absorption and emis- 
sion charge-transfer bands. Since the analytical theory 
seems to be consistent with the computer experiment, 
one needs a test against experimental data. Unfortu- 
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TABLE III: Solvation Gibbs energy (kcal/mol) and solvation 
entropy (e.u., cal K _1 mol _1 ). 



Coo 


/■ i a 
U S ,1 


r* a 
Ob, 1 


G s ,2 b 




AG S C 


AS S C 


1.0 


-236.01 


-198.30 


-183.36 


-132.68 


52.65 


65.62 


1.2 


-245.13 


-228.90 


-189.28 


-151.94 


55.85 


76.96 


1.4 


-255.10 


-262.71 


-195.75 


-173.15 


59.35 


89.56 


1.6 


-266.00 


-299.28 


-202.79 


-195.60 


63.21 


103.68 


1.8 d 


-277.88 


-337.52 


-210.43 


-219.15 


67.45 


118.37 


2.0 


-290.77 


-378.20 


-218.66 


-243.41 


72.11 


134.79 


1.8 e 


-325.1 


-2.9 


-242.1 


-2.0 


83.0 


1.1 


1.8' 


-196.0 


-1.46 


-277.88 


-1.21 


67.45 


0.26 



"Gibbs energy and solvation entropy in the initial ET state. 
6 Gibbs energy and solvation entropy in the final ET state. 
C AG S = G s ,2 — G s ,i, ASs = Ss,2 — S s ,i- 
^Parameters corresponding to water at ambient conditions. 
e Continuum calculation (DelPhi |3) with solute's vdW cavity 
^Continuum calculation (DelPhi [7]) with solute's SA cavity 

nately, experimental evidence on the solvent e ntrop ic ef- 
fects on ET reactions is very limited (see Ref. Il20l for a 
recent review). 




-4 -2 2 4 

XA 

FIG. 11: Diabatic initial (solid lines) and final (dashed lines) 
curves obtained from Eq. 1211 for parameters of DSA 1. 
Curves marked with n and m indexes refer to vibrational 
states of the initial and final states, respectively. Relative 
energies are drawn to scale, based on A s =45.15 kcal/mol, 
AG = -41.2 kcal/mol, 7kj„=400 cm" 1 . 



B. ET rate constant 

The calculations of the temperature dependent reorga- 
nization energy and equilibrium energy gap can be com- 
pared to experimental Arrhenius law measurements |4lj 
for complex 1. Transition metal-based charge-transfer 
complexes are commonly characterized by metal-ligand 
vibrational frequencies |4(j in the range u v ~ 300 — 500 
cm , substantially lower than frequencies uj v ~ 1100 — 
1500 cm -1 normally assigned to C— C skeletal vibrations 
of organic donor-acceptor complexes. Therefore, Eqs. 
JTJ, J2Jl, and (|21[1 with the full quantum-mechanical de- 
scription of vibrations and temperature-induced popu- 
lations of vibrational states should be used for the ET 
rate in complex 1. Unfortunately, our calculations pro- 
vide only the solvent component of the free energy gap. 



Its gas-phase component is unknown and the electronic 
coupling entering the Golden Rule ET rate in Eq. (JTJ 
is known with uncertainty |40j. These two parameters 
(AG g and Vn) were varied in fitting the experimental 
activation enthalpy AiJt — 9.5 kcal/mol and the experi- 
mental activation entropy AS^ /kn = —5.6 e.u. |4l|. Note 
that the experimental quantity is an effective entropy, 
including contributions due to the electronic coupling el- 
ement as well as solvation and inner-sphere vibrational 
modes 40] . The Arrhenius analysis was performed by 
the linear regression of h\(k,ET /T) vs 1/T based on the 
transition-state expression 

fc ET = ^ e -AGt(T)/ fcB T (7?) 

The rate constant was calculated based on Eqs. JTJ 
and (I21I1 . with A s and AG S varied linearly with temper- 
ature using the calculated entropies (Tables ITU and ITTTf l . 
The results of calculations are listed in Table IIVI The 
fitted electronic coupling V\i falls in the range of values 
given by electronic structure calculations [40j using the 
semi empirical INDO/s model by Zerner and co-workers 
|l2lj . The equilibrium gap obtained from the fit is ap- 
preciably more negative than AG ~ —25.4 kcal/mol es- 
timated from the redox potentials of separate donor and 
acceptor sites, based on the high spin ground state of 
the Co 2+ product (it has been argued |4£| that the less 
exothermic low spin Co 2+ product may be the relevant 
one in the experimentally observed process). Neglecting 
the vibrational excitations in the analysis (0-0 transition 
only) results in a much lower activation enthalpy and a 
substantially more negative activation entropy (second 
row in Table ITV)l . 

The relatively low frequency of metal-ligand vibrations 
in transition metal complexes results in a dense manifold 
of vibrational levels (Fig. 11 If) which are partially popu- 
lated at room temperature. The change of the vibrational 
populations with temperature may res ult i n a contribu- 
tion to the overall activation entropy |l22j . This, how- 
ever, does not happen for complex 1 when A s and AG S 
are fixed at their 298 K values. The dashed lines in Fig. 
H2l show the enthalpy and entropy of activation as a func- 
tion of the vibrational frequency at constant temperature 
and A„. Increasing the vibrational frequency makes vi- 
brational excitations less accessible, but this is seen to 
have little effect on the activation entropy and enthalpy. 

This situation changes when the temperature depen- 
dence of A s and AG S is included in the calculations of the 
Arrhenius activation parameters. In this case, the tem- 
perature dependence of the ET energy gap results in a 
change of the vibrational quantum numbers correspond- 
ing to the maximum vibronic contribution. The splitting 
of the activation barrier into the entropic and enthalpic 
contribution then becomes sensitive to the choice of lo v 
(Fig. 1121 solid lines). This sensitivity may be important 
for the interpretation of experimental data since the cor- 
rect definition of the effective vibrational frequency [Eq. 
(|15fl ] increases in importance once the temperature de- 
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FIG. 12: Enthalpy and entropy of activation of DSA complex 
1 vs ?uj v . The solid lines indicate the calculations according 
to Eq. 12111 with A s and AG S varying linearly with temper- 
ature based on the corresponding entropies from Tables UTI 
and lllll The dashed line indicate the same calculation with 
A s and AG S fixed at their 298 K values. The dash-dotted 
(Sx = 0, AS S = 0) and dotted (Sx / 0, AS S / 0) lines 
for the activation enthalpy and entropy refer to the ui v — » oo 
limit corresponding to Eq. 1131 with A„ = (no vibrational 
excitations) . 

pendence of the solvation parameters is introduced into 
the analysis of reaction rates. The classical Marcus-Hush 
equation with A„ = replaces the sum over all possible 
vibronic transitions with a single 0-0 transition. The re- 
sult is a significantly lower enthalpy and more negative 
entropy of activation (Table llV|l . 



TABLE IV: Parameters for complex 1 at T=298 K. 



Level 


Vl2 


AG 


hui v 


A-u 


AS 1 






cm -1 


kcal/mol 


cm -1 


kcal/mol 


e.u. 


kcal/mol 


Full 


0.07 a 


-41. 2 a 


400 


16.1" 


-5.5 C 


9.4 C 


A„ =0 d 


0.07 


-25.4 e 




0.0 


-25.0 d 


2.0 d 



"Obtained from fitting the experimental Arrhenius dependence. 
Sx and AS 3 are calculated within the PPSF with parameters cor- 
responding to ambient water. 

'From Ref. ^ 

Experimental values from Ref. 1411 

^Obtained by neglecting intramolecular vibrations (A„ = 0) in Eq. 
1211 . In this limiting case, the V12 value was taken from the value 
obtained from the full analysis. 

e Estimated from redox potentials of separate donor and acceptor, 
Ref. El 



VI. DISCUSSION 

The most relevant question in comparing microscopic 
solvation theories with the dielectric continuum approxi- 
mation is why the latter has allowed to describe so many 
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FIG. 13: Microscopic structure factors and their continuum 
limits. 

systems after proper parameterization of dielectric cavi- 
ties, despite drastic approximations involved. The micro- 
scopic NRFT formulation contains dielectric continuum 
as its limit, allowing us to address this question. The 
continuum limit is obtained by neglecting the spatial cor- 
relations between solvent dipoles, i.e. by neglecting the fc- 
dependence in the polarization response functions. This 
implies that polarization structure factors are replaced 
by their k = values (Fig. I13|) . This replacement is not 
a good approximation for the transverse structure factor, 
which changes quite sharply even for small k- values, but 
may be a reasonable approximation for the longitudinal 
structure factor, which is relatively flat in the range of 
fc-values significant for solvation thermodynamics. How- 
ever, for most charge configurations, even for the point 
dipole |35l ] , the contribution of transverse polarization to 
the solvation free energy is relatively small [HI] (~ 10% in 
our calculations for complex 1 in water). Therefore, the 
inaccurate continuum approximation for the transverse 
structure factor does not significantly affect the results 
of calculations. 

The continuum estimates for the polarization structure 
factors result in the following inequalities between the 
continuum and microscopic longitudinal and transverse 
components of the reorganization energy 

A L,co„t < A ^ A T,co„t > A T ( 7g ) 

The sharp change of the transverse structure factor at 
small fc-values is responsible for a substantial overesti- 
mate of the transverse component of solvation by contin- 
uum models [35, 36] . This overestimate manifests itself in 
solvation dynamics. The transverse polarization dynam- 
ics is much slower than the longitudinal polarization dy- 
namics |27j . Therefore, continuum models predict bipha- 
sic solvation dynamics with an appreciable slow compo- 
nent due to transverse polarization relaxation. This slow 
component is not o bser ved in computer simulations of 
solvation dynamics 123] and it does not show up in the 
microscopic calculations reported in Ref. |3!| 

The relatively flat form of the longitudinal structure 
factors at low fc-values does not mean that replacing 
S L ' T (k) by S L,T (0) gives accurate numbers for the sol- 
vation free energy and/or the reorganization energy. A 
moderate increase of S L (k) in the range of wavevectors 
contributing to the fc-integral substantially affects the 
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FIG. 14: Reorganization energies in polar solvents vs the 
Pekar factor [Eq. I0HJ] (a) and the Lippert-Mataga polar- 
ity parameter [Eq. JZHJ] (b). The open points indicate the 
DelPhi calculations with vdW (squares) and SAS (circles) 
molecular surfaces. The closed points (triangles) refer to 
calculations with the NRFT. Numbers on the plot indicate: 
chloroform (1), tetrahydrofuran (2), methylacetate (3), N,N- 
dimethylformamide (4), acetone (5), acetonitrile (6), water 
(7). 



calculated values of solvation free energies (cf. columns 5 
and 6 in Table HJ). Moreover, the gap between the micro- 
scopic and continuum values changes with the solvent di- 
electric parameters (see, e.g., Fig. EU- This observation 
practically means that there is fundamentally no unique 
scheme for defining the dielectric cavity applicable to all 
solvent polarities. 

The dominance of longitudinal polarization fluctua- 
tions in solvation thermodynamics is also responsible for 
experimentally observed linear trends o f the reorgani- 
zation energy with the Pekar factor p^l Il24l Il25j | [Eq. 
I|48p] . Even at the continuum level, the polarization re- 
sponse function for a solute of complex shape is not repre- 
sented by the Pekar factor appearing in the longitudinal 
projection of the solvent response function |l2fij . How- 
ever, large separation of charges is responsible for the 
predominantly longitudinal response of the solvent, and 
continuum reorganization energies calculated for complex 
1 in polar solvents correlate well with the Pekar factor 
(Fig. 114b ). If fact, an equally good correlation is seen in 
respect to the Lippert-Mataga polarity parameter com- 
monly used for solvation of dipoles (Fig. 114b): 



/•LM _ 
J n 



e.,-1 



1 



2c, 



1 2e, 



1 



(79) 



The use of a particular parameter does not therefore tell 
much about the nature of the solute charge distribution 
and, obviously, reflects a linear relation between cq and 
f^ M for common solvents. 

The results of the current microscopic calculations are 
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FIG. 15: Reorganization entropies from continuum calcula- 
tions (DelPhi with the vdW surface, closed circles) and by 
the NRFT (open squares). Points refer to the same solvents 
as in Fig. 1141 



shown by triangles in Fig.^J These numbers do not ex- 
hibit a linear dependence, although the extent of scatter 
is not uncommon for ET experiment. The comparison of 
the continuum and microscopic dependence on the sol- 
vent polarity does not permit a clear distinction between 
the two formulations. Where the distinction becomes 
clear is for the reorganization entropy in strongly polar 
solvents. Figure 1151 shows reorganization entropies S\ 
calculated in continuum (DelPhi pfl Poisson-Bolzmann 
solver) and microscopic (NRFT) theories. The contin- 
uum calculation reflects the temperature variation of the 
Pekar factor cq: 



-(dc /dT)i 



\de s /dT) P (80) 



In low-polarity solvents, Co is mostly influenced by the 
static dielectric constant, which has a negative temper- 
ature derivative. The continuum reorganization entropy 
(closed circles in Fig. lFo|) is positive and is close to the 
microscopic result (open squares in Fig. 115ft . The contin- 
uum estimate of the temperature variation of A s in low- 
polarity solvents thus gives a semi-quantitative account 
of the experimental observations 51] . In strongly polar 
solvents, the temperature derivative of Cq is mostly influ- 
enced by the high-frequency dielectric constant, and con- 
tinuum S\ is nagative. In this case, the predictions of the 
continuum model significantly depart from both the mi- 
croscopic calculations and many experimental measure- 
ments [12J, |4jJ |4J, |4jJ |4y, |47], |50j , showing positive reor- 
ganization entropies. 

The microscopic calculations presented here show a rel- 
atively weak dependence of the reorganization energy on 
the solvent high-frequency dielectric constant, in qual- 
itative accord with available computer simulation data 
[ftL ll8L Il9| . Testing this theoretical prediction experi- 
mentally may become problematic because of the narrow 
range of values available for common polar solvents. 
We note, however, that the problem of the weak depen- 
dence of the reorganization energy on is related to 
the problem of correct sign of the reorganization entropy. 
The strong dependence of the continuum reorganization 
energy on f M is one of major factors shifting the con- 
tinuum reorganization entropy to the range of positive 
values. 
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The calculations of the quadrupolar component of 
the solvent reorganization energy presented here con- 
firm the conclusion previously reached for Stokes shifts in 
coumarin-153 optical dye jl27| : quadrupolar solvation is 
insignificant in most commonly used polar solvents, and 
the dipolar approximation for the solvent charge distri- 
bution is sufficient for most practical calculations. 



Acknowledgments 

D.V.M. thanks the Donors of The Petroleum Research 
Fund, administered by the American Chemical Society 
(39539-AC6), for support of this research. M.D.N, was 
supported by DE-AC02-98CH10886 at Brookhaven Na- 
tional Laboratory. The authors are grateful to Prof. 
G. A. Voth for sharing the structural data on the 
polypeptide-linked donor-acceptor complex. This is pub- 
lication #596 from the ASU Photosynthesis Center. 



APPENDIX A: SIMULATION AND ANALYSIS. 

The MC simulations of dipolar-polarizable hard sphere 
solvents shown in Fig. [7| were done as described in Ref. 
Il9l Simulations of 6 x 10 5 cycles long were run for 1372 
polarizable molecules with periodic boundary conditions 
and the reaction field cutoff of dipole-dipolc interactions. 
The MD simulations were carried out with the force field 
of 3-site acetonitrile (ACN3) by Edwards, Madden, and 
McDonald [lO^ and the 3-site model of water (TIP3P) by 
Jorgensen et al. |lOlj ( Table fvjl. The site-site interaction 
potential is given by the sum of the Lennard-Jones (LJ) 
and Coulomb interaction potentials: 



<Jgf3 

r a /3 



Tap 



gag/3 

r a /3 



(Al) 



where the LJ parameters are taken according to the 
Lorentz-Bertholet rules: e a/ 3 = ^Je a £p and a a p = (a a + 
cp)/2. All simulations were don e with the DL_POLY 
molecular dynamics package |l28j . We run two sets of 
MD simulations in the temperature range from 288 K to 
308 K with a 5 K step. The timestep in each simulation 
is 5 fs. All MD simulation are 20 ns long. 



TABLE V: Force field parameters for use in the simulation. 
Atomic interaction site a a /A e a x 1Q A / (kcal/mol) q a j e 



TIP3P water" 








O 


3.15 


152.10 


-0.834 


Acetonitrile 6 








Me 


3.6 


379.55 


0.269 


C 


3.4 


99.36 


0.129 


N 


3.3 


99.36 


-0.398 



"roH=0.9572 A, Z HOH=104.52° 
i, r M eC=l-46 A, r CN =1.17 A 

We used the Nose-Hoover thermostat |129| | for the 
ACN3 simulations with the relaxation parameter of 0.5 
fs. This value ensures good stabilization of the total sys- 
tem energy. The energy drift for ACN3 is only about 
0.1%. The simulation box was constructed to include 
256 ACN3 molecules in a cube with the side length 
L = 28.2025 A at T=298 K to reproduce the experi- 
mental mass density of acetonitrile, pAf=0.777 g/cm 3 . 
The side length is adjusted at each temperature to ac- 
count for temperature expansion with the experimental 
volume expansion coefficient a p — 1.38 x 10~ 3 K . 

In simulations of TIP3P water, 256 molecules reside in 
a cube with the side length of L = 19.7744 A at 298 K. 
The system is coupled to the Berendsen |13C| thermostat 
with the relaxation time of 0.1 fs. The drift in total en- 
ergy of about 0.1 % is observed. The liquid mass density 
Pm = 0.9896 g/cm 3 and the volume expa nsion coefficient 
a p = 2.96 x 10~ 3 K" 1 are taken from Ref.lllcL The latter- 
value is close to the experimental expansion coefficient of 
ambient water, a p = 2.6 x 10~ 3 K _1 . 

The cutoff for short-range LJ interaction is 13 A for 
ACN3 and 9 A for TIP3P. For long-range Coulomb inter- 
actions, Ewald summation from DL_POLY |l3lj | is used 
for ACN3 and smoothed particle mesh (SPME) |l3^ 
Ewald is adopted for TIP3P. Ewald summation parame- 
ters are the convergence parameter a and the maximum 
wavenumber k™y X z . The parameter sets a = 0.24 A -1 , 

= 7 A-1 >' and a = - 35 A_1 > = 8 A" 1 were 

used for ACN3 and TIP3P respectively. 

The structure factors have been calculated as the vari- 
ance of longitudinal and transverse projections of the It- 
space solvent polarization 



JV 



M(fe) = (l/m)^m ie - 



(A2) 



where m; = J2 a 9 ar ? ^ s a dipole moment of the ith 
molecule and the sum runs over the N molecules in the 
simulation box. The static dielectric constant is given in 
terms of the k = variance as follows Il33ll 



l + 3y(M(0) 2 )/N, 



(A3) 



where y = (47r / '9) pm 2 / k B T . 
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